# 已知LV模型的参数，求变化周期等性质。
# 首先载入数值计算模块 numpy, 符号计算模块 sympy, 画图模块 pylab, 和 科学计算的积分模块 scipy.integrate. 

import numpy as np
import sympy as sp
import pylab as plt
from scipy.integrate import odeint

# 定义两个符号变量 x,y, 求解方程组，使用符号模块的 solve 函数计算平衡点。

x=sp.symbols('x')
y=sp.symbols('y')
eq=[0.2*x-0.005*x*y,-0.5*y+0.01*x*y]
s1=sp.solve(eq,[x,y])
print(s1)

# 在 x,y 平面的正方形区域 $[0,100]\times [0,100]$ 上画出格点，间隔距离为10, 共121个点。使用 quiver 函数画出方向场。结果见图1. 

x=np.linspace(0,100,11)
x,y=np.meshgrid(x,x)
u=0.2*x-0.005*x*y; v=-0.5*y+0.01*x*y
plt.quiver(x,y,u,v)
plt.xlabel('$x$')
plt.ylabel('$y$',rotation=0)

# 将常微分方程的右端创建为一个函数，设自变量取值区间为 $[0,100]$, 自变量按间隔 $\Delta t =0.1$ 设置分点 $t_i, 0\le i\le 1000$. 使用 odeint 函数求解常微分方程，数值解 $x(t_i), y(t_i)$ 保存在变量 s 中。

def func(f,t):
    x, y=f
    return [0.2*x-0.005*x*y,-0.5*y+0.01*x*y]
t=np.linspace(0,100,1000)
s=odeint(func, [70,40], t)

# 求出被捕食者的总体数量 $x(t)$ 的最大值和最小值，捕食者的总体数量 $y(t)$ 的最大值和最小值。

x1=max(s[:,0])
x2=min(s[:,0])

print('The maximal value of x is ', x1)
print('The minimal value of x is ', x2)

y1=max(s[:,1])
y2=min(s[:,1])

print('The maximal value of y is ', y1)
print('The minimal value of y is ', y2)

# 在前面用 quiver 函数画的方向场中画出解函数的轨线 $(x(t),y(t))$. 结果见图1. 

plt.plot(s[:,0], s[:,1])

# 另起一个图，以时间 $t$ 为横坐标，画出被捕食者的总体数量 $x(t)$ 的函数图像，与捕食者的总体数量 $y(t)$ 的函数图像。结果见图2. 

plt.figure()
plt.plot(t, s[:,0], 'r-', label='$x(t)$')
plt.plot(t, s[:,1], 'b--', label='$y(t)$')
plt.legend()
plt.xlabel('$t$')

